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Recent  climate  changes  have  led  to  an  increase  in  the  exposure  of  archaeological  remains  in  frozen 
environments  due  to  the  melting  of  glaciers  and  ice  patches,  and  the  thawing  of  permafrost.  In  some 
cases,  the  discovery  of  glacial  archaeological  findings  has  occurred  due  to  chance.  In  order  to  avoid  the 
risk  of  losing  exceptional,  often  organic,  cultural  remains  due  to  decomposition,  systematic  and  pre¬ 
dictive  methods  should  be  employed  to  locate  areas  of  high  glacial  archaeological  potential.  Here,  we 
merged  archaeological  and  glaciological  methods  to  create  a  new  type  of  archaeological  prediction 
model  in  the  field  of  glacial  archaeology.  Locational  analysis  and  glaciological  modelling  were  used  to 
highlight  current  and  future  areas  of  archaeological  potential  in  the  Pennine  Alps,  located  between 
Switzerland  and  Italy.  Future  glacier  area  was  calculated  in  10  year  increments  until  2100.  By  2090,  93%  of 
glacier  area  is  expected  to  have  disappeared.  The  results  from  the  final  model,  GlaciArch,  provide  new 
insights  into  future  glacial  archaeological  prospection  in  the  Pennine  Alps  by  narrowing  down  a  study 
region  of  4500  km2  into  several  manageable  square  kilometre  sites. 

©  2014  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Due  to  the  alternation  of  various  warm  and  cold  periods,  glacier 
extents  and  ice  volume  storage  have  fluctuated  in  the  entire  Eu¬ 
ropean  Alps  during  the  Holocene  (10.5  ka  to  present).  Compared  to 
the  Last  Glacial  Maximum  (LGM)  (19—20  ka  BP)  and  latest  Pleis¬ 
tocene,  when  large  piedmont  lobes  of  vast  valley  glaciers  reached 
the  Alpine  foreland  (Clark  et  al.,  2009;  Ivy-Ochs  et  al„  2008),  glacier 
changes  have  been  rather  minor  during  the  Holocene.  The  glaci- 
erized  area  varied  between  the  stage  of  the  Little  Ice  Age  (LIA) 
maximum,  around  1850,  and  a  minimum  which  was  significantly 
smaller  than  the  present  day  extents  (Grosjean  et  al.,  2007; 
Holzhauser,  2007;  Joerin  et  al„  2006,  2008). 

Glacier-climate  interactions  have  affected  humans  for  millennia. 
In  the  European  Alps,  glacier  fluctuations  directly  influenced  hu¬ 
man  interaction  with  Alpine  areas  (Benedict  and  Olson,  1978; 
Wiegandt  and  Lugon,  2008).  For  example,  as  glaciers  receded  af¬ 
ter  the  LGM,  humans  took  advantage  of  the  newly  ice-free  Alpine 
biome  which  offered  plenty  of  food  and  resources  during  the 
Paleolithic  period  (Pacher,  2003;  Tagliacozzo  and  Fiore,  2000).  The 
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present  atmospheric  warming  has  caused  shrinkage  of  glaciers  and 
ice  caps  all  over  the  world  (IPCC,  2013).  In  consequence,  melting  ice 
and  snow  has  uncovered  archaeological  remains  in  Arctic  and 
Alpine  environments  (Andrews  et  al„  2012;  Beattie  et  al.,  2000; 
Callanan,  2012,  2013;  Dixon  et  al„  2005;  Farbregd,  1972;  Farnell 
et  al„  2004;  Hafner,  2012;  Hare  et  al.,  2004,  2012;  Lee,  2012; 
Rogers  et  al.,  2014;  VanderHoek  et  al„  2007)  which  further  attests 
to  the  use  of  frozen  regions  on  a  global  scale.  These  artefacts  which 
have  melted  out  of  ice  patches  and  glaciers,  and  thawed  out  of 
permafrost,  have  created  a  new  sub-discipline  of  archaeology: 
glacial  archaeology.  “Glacial  archaeology”  has  also  been  referred  to 
as  ice  patch  archaeology  (c.f.  Andrews  and  MacKay,  2012;  Reckin, 
2013)  and  frozen  archaeology  (Molyneaux  and  Reay,  2010). 
Perhaps  one  of  the  most  famous  examples  of  a  glacial  archaeolog¬ 
ical  find  is  that  of  Otzi  the  Tyrolean  Iceman  who  was  accidentally 
discovered  by  hikers  in  1991  on  the  Italian/Austrian  border,  pro¬ 
truding  from  an  ice  patch  (Prinoth-Fornwagner  and  Niklaus,  1994; 
Seidler  et  al„  1992).  The  uniqueness  of  Otzi  and  other  glacial 
archaeological  discoveries  is  that  they  have  often  been  preserved  by 
ice  for  thousands  of  years,  thus  protecting  them  and  providing 
scientists  with  unparalleled  information  about  past  cultures  and 
climates  (Dixon  et  al.,  2005;  Reckin,  2013).  There  is  urgency  to 
collect  these  delicate,  often  organic,  glacial  archaeological  remains 
before,  or  soon  after,  they  melt  out  of  the  ice  and  become  destroyed 
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by  decomposition  (Andrews  and  MacKay,  2012;  Dixon  et  al.,  2005; 
Molyneaux  and  Reay,  2010).  As  melting  in  high  altitudes  and  lati¬ 
tudes  is  not  anticipated  to  halt  in  the  near  future  (c.f.  Radic  et  al., 
2014),  more  glacial  archaeological  finds  can  be  expected  and 
there  is  a  need  to  further  develop  predictive  methods  in  this 
research  domain. 

In  this  paper,  archaeological  and  glaciological  methods  are 
merged  together  to  create  a  new  type  of  predictive  model  to 
determine  areas  of  glacial  archaeological  potential.  The  results  will 
be  used  as  a  decision  support  tool  for  future  prospection  of 
archaeological  findings  in  high  mountain  environments.  Our 
approach,  referred  to  as  “GlaciArch”  in  the  following,  is  based  on 
current  ice  thickness  distribution,  future  evolution  of  glacierized 
areas,  and  topographic  characteristics  of  the  terrain  which  could 
have  influenced  past  human  accessibility.  First,  currently  glaci¬ 
erized  or  recently  deglacierized  high  altitude  mountain  passes 
located  on  the  border  of  Switzerland  and  Italy  are  selected  to  be 
used  as  sites  on  which  to  perform  the  locational  analysis.  Next,  least 
cost  paths  (LCPs)  are  calculated  between  valleys  and  respective 
passes.  Then,  locational  analysis  is  used  to  determine  areas  of 
glacial  archaeological  potential  based  on  the  physical  characteris¬ 
tics  of  the  terrain.  After,  the  future  evolution  of  glaciers  is  modelled 
for  the  Pennine  Alps  using  a  glacier  evolution  model  (Huss  et  al., 
2010a).  Finally,  the  results  of  glacier  modelling  are  combined 
with  the  results  of  locational  analysis  to  create  GlaciArch,  a  pre¬ 
dictive  model  which  ultimately  defines  regions  of  highest 


archaeological  interest  for  now  and  the  future.  This  paper  high¬ 
lights  how  the  intersection  of  glaciological  and  archaeological 
methods  provides  a  new  approach  for  looking  at  glacial  archaeo¬ 
logical  prospection. 

2.  Study  area  and  data 

2.3.  Study  area 

The  Pennine  Alps  (centered  at  approximately  45°57'N,  7°32'E) 
are  located  between  the  canton  of  Valais,  Switzerland,  and  the 
provinces  of  Aosta  and  Piedmont,  Italy  (Fig.  1).  The  whole  region  is 
of  particular  glacial  archaeological  interest  due  to  its  large  glaci¬ 
erized  area  and  rich  cultural  heritage.  The  Pennine  Alps  cover 
approximately  4500  km2  and  reach  altitudes  above  4000  m  a.s.l. 
The  main  valleys  to  the  north,  south,  and  east  of  the  Pennine  Alps, 
the  Rhone  valley  (Switzerland),  and  the  Aosta  and  Antigorio  valleys 
(Italy)  respectively,  are  scattered  with  archaeological  remains 
dating  from  Mesolithic  (9.5  ka  to5.5  ka  BC)  to  historic  times  (Curdy, 
2007;  Radmilli,  1963).  Although  most  travellers  reached  these 
valleys  from  lower  altitudes,  each  valley  could  also  be  reached  by 
crossing  the  Pennine  Alps  between  them.  This  relatively  short 
distance  was  often  traversed  for  commercial  purposes.  Archaeo¬ 
logical  remains  collected  on  the  way  to,  and  on  top  of,  mountain 
passes  between  Switzerland  and  Italy  demonstrate  the  use  of  these 
passes  as  trade  and  travel  routes  for  thousands  of  years  (Bezinge 
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Fig.  1.  Overview  of  study  area.  Glacierized  areas  are  shaded  in  dark  grey. 
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and  Curdy,  1994,  1995;  Coolidge,  1912;  Curdy,  2007;  Curdy  et  al., 
2003;  Harriss,  1970,  1971;  Lehner  and  Julen,  1991;  Rogers  et  al., 

2014).  In  fact,  early  Neolithic  culture  seems  to  have  spread  to 
Valais  via  the  high  altitude  passes  of  the  Pennine  Alps  from  the 
south,  possibly  due  to  the  grazing  of  small  herds  in  the  high  pas¬ 
tures  in  summer  (Curdy  et  al.,  2003;  Curdy,  2007).  Throughout 
Prehistory  and  the  Roman  Period,  there  were  indications  of  a  strong 
cultural  relationship  between  Aosta  and  the  Rhone  valley;  later,  the 
Aosta  and  Upper  Rhone  valleys  were  integrated  as  the  unique 
ecclesiastical  province  of  Tarentasia  for  several  centuries  (Harriss, 
1970;  Curdy,  2010).  The  relatively  few  archaeological  remains 
found  at  high  altitudes  in  this  region  should  not  be  considered  to  be 
a  direct  result  of  the  use  of  these  high  altitude  passes.  In  the  past, 
contrary  to  current  beliefs,  high  altitude  regions  were  used  more 
often  than  assumed,  and  proved  to  be  more  hospitable  than  they 
seem  to  modern  day  people  (Aldenderfer,  2006;  Reckin,  2013; 
Walsh  et  al.,  2006). 

2.2.  Data 

The  high  altitude  pass  names  and  locations  used  in  the  first  step 
of  the  locational  analysis  are  derived  from  the  25  m  resolution 
SwissNames  database  provided  by  the  Swiss  Federal  Office  of 
Topography  (swisstopo)  (Federal  Office  of  Topography  (2014)) 
which  contains  all  names  given  on  the  1:25,000  national  topo¬ 
graphic  maps.  The  1973  Swiss  glacier  inventory  (Muller  et  al.,  1976) 
was  used  to  determine  glacierized  or  recently  deglacierized  passes. 
The  Digital  Elevation  Model  (DEM)  used  in  both  the  Least  Cost  Path 
Analysis  (LCPA)  and  slope  calculation  was  the  global  30  m  resolu¬ 
tion  Advanced  Spaceborne  Thermal  Emission  and  Reflection  Radi¬ 
ometer  Global  Digital  Elevation  Model  (ASTER  GDEM)  (version  2) 
(NASA,  2012).  Although  more  accurate  DEMs  exist,  it  was  not 
possible  to  obtain  a  consistent  DEM  for  each  side  of  the  Pennine 
Alps.  The  ASTER  GDEM  provides  a  consistent  and  sufficient  data 
accuracy  for  these  regional  scale  calculations  in  this  study. 

For  the  Swiss  glaciers,  the  new  Swiss  Glacier  Inventory  SGI2010 
(Fischer  et  al.,  in  press)  was  used.  This  layer  was  created  by  manual 
digitization  from  high  resolution  (50  cm)  aerial  orthoimagery  ac¬ 
quired  between  2008  and  2011.  For  the  Italian  glaciers,  outlines 
based  on  satellite  imagery  of  2003  were  used  (Paul  et  al.,  2011). 
Surface  topography  for  each  glacier  was  extracted  by  intersecting 
glacier  outlines  with  terrain  elevation  data.  Glacier  ice  thickness 
distribution  and  bedrock  topography  were  calculated  based  on  a 
flux-gate  approach  and  using  the  principles  of  ice  flow  dynamics 
(Huss  and  Farinotti,  2012).  Information  on  surface  mass  balance  of  a 
large  sample  of  glaciers  in  the  Pennine  Alps  over  the  last  decades  is 
available  from  a  combination  of  direct  field  observations,  geodetic 
ice  volume  changes  and  distributed  modelling  (Huss,  2012).  Sce¬ 
narios  for  the  future  evolution  of  climatological  variables  are  ob¬ 
tained  from  Regional  Climate  Models  (RCM)  from  the  CH2014 
project  (CH2014-Impacts,  2014). 

3.  Background 

In  the  Swiss  Alps,  glacial  archaeological  finds  have  been  located 
at  four  locations:  two  sites  on  the  border  of  the  cantons  of  Valais 
and  Bern  in  the  Bernese  Alps,  the  Lotschenpass  (Bellwald,  1992; 
Meyer,  1992)  and  Schnidejoch  pass  (Hafner,  2012);  one  site  in 
eastern  Switzerland,  the  Porchabella  glacier  (Rageth,  1995);  and  in 
Valais,  located  in  the  Pennine  Alps,  near  the  Theodulpass  (Lehner 
and  Julen,  1991;  Meyer,  1992).  The  oldest  and  most  notable  finds 
were  discovered  between  2004  and  2011  at  the  Schnidejoch  pass 
from  a  melting  ice  patch  which  was  formerly  attached  to  the 
Chilchli  glacier  on  the  north  side  of  the  pass.  The  ages  of  the  finds 
range  from  the  Neolithic,  Early  Bronze  Age,  Iron  Age,  Roman,  and 


Medieval  periods  making  this  one  of  the  most  prolific  glacial 
archaeological  sites  in  the  Alps  (Hafner,  2012).  The  abundance  of 
finds  can  be  attributed  to  the  location  of  the  ice  patch  which  is  in  a 
small  depression  facing  northeast  where  ice  has  accumulated  over 
centuries.  The  Pennine  Alps'  most  prolific  glacial  archaeological  site 
to  date  has  been  that  of  the  Theodulpass.  The  finds,  which  include 
skeletal  remains,  leather  clothing  and  shoe  soles,  weapons,  and 
coins  from  the  “Mercenary  of  Theodul",  date  back  to  the  16th 
century  (Lehner  and  Julen,  1991;  Meyer,  1992).  These  items  were 
found  between  1985  and  1990  along  the  margins  of  the  Oberer 
Theodul  glacier  on  the  Swiss  side  of  the  border.  It  is  believed  that 
the  Mercenary  fell  into  a  crevasse  and  was  preserved  for  hundreds 
of  years  until  glacial  dynamics  and  melting  eventually  released  him 
and  his  belongings. 

3.1.  Archaeological  predictive  modelling 

In  archaeology,  the  use  of  predictive  modelling  began  in  the 
1980's  and  has  since  grown  into  a  wide  research  field,  mostly  due  to 
an  increase  in  the  accessibility  to  Geographic  Information  Systems 
(GIS)  software  and  the  ever-improving  spatial  resolution  of  data 
(c.f.  Ebert,  2004;  Kvamme,  1999;  McCoy  and  Ladefoged,  2009).  In 
most  cases,  archaeological  predictive  modelling  is  used  to  forecast 
the  location  of  archaeological  sites  based  on  the  presence  or 
absence  of  defined  criteria,  in  order  to  allocate  information  about 
known  patterns  onto  unknown  places  (Conolly  and  Lake,  2006; 
Warren  and  Asch,  2003;  Wheatley  and  Gillings,  2002).  Inputs  to 
predictive  models  usually  include  sampled  sites  and  have  the  main 
goal  of  finding  new  sites  which  were  used  for  human  occupation 
(Carleton  et  al.,  2012;  Carrer,  2013;  Graves,  2011 ;  Kohler  and  Parker, 
1986).  Put  simplistically,  predictive  methods  have  been  used  to 
determine  the  level  of  archaeological  potential  in  a  region  and 
provide  decision-makers  with  a  tool  to  justify  why  certain  areas  are 
more  archaeologically  interesting  than  others  (McCoy  and 
Ladefoged,  2009). 

In  glacial  archaeology  specifically,  predictive  methods  have  been 
used  in  relatively  few  instances  but  show  promising  results 
(Andrews  et  al.,  2012;  Dixon  et  al.,  2005).  For  example,  in  Alaska, 
Dixon  et  al.  (2005)  were  the  first  to  use  predictive  modelling  for 
glacial  archaeological  purposes  by  using  a  weighted  combination  of 
cultural,  biological,  and  geological  input  layers  to  successfully 
determine  areas  of  high  glacial  archaeological  potential.  Similarly, 
Andrews  et  al.  (2012 )  used  remotely  sensed  data  and  other  weighted 
input  layers  to  determine  areas  of  glacial  archaeological  potential  in 
northern  Canada.  Both  previous  studies  were  conducted  at  various 
ice  patch  sites.  This  study,  which  focuses  on  the  potential  of  locating 
glacial  archaeological  remains  on  or  near  glaciers,  differs  from  pre¬ 
vious  ones  based  on  the  distinctive  environments.  Ice  characteristics 
and  glacier  dynamics  can  strongly  affect  the  potential  of  locating 
glacial  archaeological  remains,  and  thus  should  be  researched  ac¬ 
cording  to  local  conditions.  For  example,  glaciers  composed  of  thick 
ice  or  located  on  steep  slopes  move  relatively  quickly  and  would 
destroy  anything  entrained  within  it  in  a  matter  of  a  few  hundred 
years  based  on  the  principles  of  ice  dynamics  (Benn  and  Evans,  2010; 
Dixon  et  al.,  2005;  Hafner,  2012).  The  margins  of  slower-moving 
glaciers  could  prove  to  be  a  better  environment  for  finding  glacial 
archaeological  remains,  with  the  ideal  environment  being  sur¬ 
rounded  by  ice  with  little  or  no  movement  such  as  ice  patches  (like 
the  Schnidejoch  site)  or  slow-moving  small  glaciers  located  on 
relatively  flat  terrains  (like  the  Theodul  site). 

3.2.  Least  cost  path  analysis 

LCPA  is  one  type  of  archaeological  prediction  method  used  in 
GIS  to  calculate  the  “optimal”  path  across  a  landscape  based  on  one 
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or  more  predefined  input  criteria  (Anderson  and  Gillam,  2000;  Bell 
and  Lock,  2000;  Egeland  et  al.,  2010;  Gorenflo  and  Gale,  1990; 
Howey,  2007;  Madly  and  Rakos,  1996;  Rogers  et  al.,  2014; 
Verhagen  and  Jeneson,  2012).  It  allows  archaeologists  to  gain  a 
better  understanding  about  movement  patterns  in  prehistoric  or 
historic  terrains  (Llobera  et  al„  2011;  Murrieta-Flores,  2010,  2012; 
White  and  Surface-Evans,  2012).  It  is  based  on  the  principal  that 
humans  will  take  the  easiest  path  from  one  location  to  another  if 
there  are  no  other  social  or  cultural  forces  directing  them  other¬ 
wise.  The  concept  is  not  unique  to  archaeology  and  was  developed 
firstly  in  psychology  and  has  since  been  used  in  various  research 
fields  (Zipf,  1949). 

3.3.  Locational  analysis 

Locational  analysis,  also  referred  to  as  archaeological  location 
modelling  (ALM)  or  site  predictive  modelling,  is  a  predictive 
method  which  calculates  archaeological  potential  based  on  multi¬ 
ple  weighted  inputs,  often  including  known  archaeological  site 
locations  (Andrews  et  al.,  2012;  Carleton  et  al.,  2012;  Carrer,  2013; 
Dixon  et  al„  2005;  Egeland  et  al.,  2010).  Like  the  term  predictive 
modelling,  locational  analysis  has  various  meanings  and  a  defini¬ 
tion  which  has  developed  over  time  (Kvamme,  1999;  McCoy  and 
Ladefoged,  2009).  Therefore,  we  define  locational  analysis  using 
McCoy  and  Ladefoged’s  (2009)  description  whereby  the  potential  of 
undiscovered  sites  will  be  determined  by  calculating  zones  of 
future  prospection  without  spatially  analysing  known  site  loca¬ 
tions.  This  method  is  often  used  in  vast  areas  from  which  archae¬ 
ological  remains  are  sparse,  possibly  due  to  a  lack  of  prospection. 

3.4.  Glacier  retreat  modelling 

The  atmospheric  warming  observed  during  the  last  decades  has 
caused  a  considerable  reduction  in  area  and  mass  of  glaciers  and  ice 
caps  all  around  the  globe  (Zemp  et  al.,  2009).  As  an  immediate 
response,  changes  in  the  climatic  forcing  acting  on  glaciers  lead  to 
changes  in  the  surface  mass  balance,  that  is,  to  changes  in  the 
quantity  of  snow  and  ice  added  to  or  melted  from  the  glacier.  The 
effects  of  the  observed  glacier  changes  are  numerous  and  apply  to  a 
broad  range  of  spatio-temporal  scales,  from  global  sea  level  rise 
(Gardner  et  al.,  2013)  to  regional  impacts  on  runoff  in  major  river 
catchments  (Kaser  et  al.,  2010),  to  local  consequences  for  landscape 
evolution,  hydropower  production,  natural  hazards  and  tourism 
(Cannone  et  al„  2008;  Farinotti  et  al.,  2012;  Fischer  et  al.,  2011). 
Consequently,  the  projection  of  future  glacier  evolution  has  gained 
increasing  attention  in  glaciology. 

In  recent  years,  various  glacier  modelling  approaches  have  been 
developed.  For  individual  glaciers,  they  range  from  simple  2D 
flowline  models  (Oerlemans  et  al.,  1998;  Van  de  Wal  and  Wild, 
2001)  to  complex  3D  coupled  mass-balance  ice-flow  models 
(Jouvet  et  al.,  2011;  Schneeberger  et  al.,  2003).  To  calculate  glacier 
response  at  the  regional  scale,  simpler  approaches  neglecting 
transient  changes  and  the  effect  of  ice  dynamics  were  used  (Paul 
et  al.,  2007;  Schaefli  et  al.,  2005).  More  advanced  models  driven 
by  distributed  surface  mass  balance  input  and  employing  a 
parameterization  of  glacier  ice  flow  have  also  been  applied  both  at 
the  single-glacier  scale  and  to  the  entire  European  Alps  (Huss, 
2012;  Huss  et  al.,  2010a;  Salzmann  et  al.,  2012). 

4.  Methods 

GlaciArch  is  composed  of  various  steps  which  culminate  in  the 
creation  of  a  predictive  model  which  gauges  the  glacial  archaeo¬ 
logical  potential  of  the  Pennine  Alps.  The  respective  steps  are 
presented  in  the  next  sections. 


4.1.  High  altitude  pass  selection  and  LCPs 

The  first  step  was  to  determine  which  high  altitude  passes  on 
the  Swiss/Italian  border  were  glacierized  in  1973.  This  coincides 
with  the  notion  that  recently  deglacierized  passes  have  a  higher 
glacial  archaeological  potential  and  should  be  prospected  first 
based  on  the  fragility  of  glacial  archaeological  remains  and  arte¬ 
facts.  The  1973  glacier  inventory  (Muller  et  al.,  1976)  was  chosen 
because  it  covers  an  ideal  time  frame  as  archaeological  items 
exposed  over  the  last  40  years  might  have  a  better  chance  to  persist 
compared  to  160  years  if  using  the  glacier  inventory  from  orthoi¬ 
magery  prior  to  the  SGI2010  (Fischer  et  al.,  2014).  The  passes  were 
chosen  using  the  selection  tools  in  ArcGIS  10.1.  First,  the  select  by 
attributes  tool  was  used  to  query  all  mountain  passes  which  could 
be  crossed  by  either  foot  or  by  road  from  the  SwissNames  database 
for  the  canton  of  Valais,  resulting  in  the  selection  of  670  records. 
Next,  the  select  by  location  tool  was  used  to  query  those  passes 
which  were  glacierized  in  1973  from  the  currently  selected  records, 
to  highlight  the  ones  that  were  recently  deglacierized,  or  still 
currently  glacierized,  leaving  111  passes.  From  those  111  passes,  the 
ones  on  and  near  the  border  between  Switzerland  and  Italy  were 
selected  as  it  is  known  that  some  of  the  high  altitude  passes  in  this 
study  region  have  been  used  for  thousands  of  years.  This  resulted  in 
the  final  extraction  of  19  border  passes  which  were  glacierized  in 
1973  from  which  to  calculate  LCPs  (Table  1,  Fig.  2a). 

LCPs  were  calculated  using  the  method  described  by  Rogers 
et  al.  (2014)  from  each  of  the  19  passes  to  their  nearest  respective 
main  valleys  (Rhone,  Aosta,  or  Antigorio)  (Fig.  2b).  This  method 
used  Tobler's  (1993)  hiking  algorithm  to  calculate  walking  times 
based  on  slope  and  prehistoric  landcover  as  inputs  (Bell  and  Lock, 
2000;  Gorenflo  and  Gale,  1990;  Rogers  et  al.,  2014;  Tobler,  1993; 
Verhagen  and  Jeneson,  2012;  Whitley  and  Hicks,  2003).  The  path 
distance  and  cost  distance  tools  in  ArcGIS,  which  calculate  the 
accumulative  cost  across  the  terrain  from  a  starting  location  and 
the  shortest  path  from  a  destination  back  to  starting  location, 
respectively,  were  used  for  the  calculations. 

4.2.  Locational  analysis 

In  this  part  of  the  analysis,  multiple  criteria  were  used  to  locate 
areas  of  high  archaeological  potential  by  analysing  where  people 
were  able  to  travel  based  on  the  topographic  characteristics  of  the 
terrain  by  measuring  the  distance  from  LCPs  and  the  slope  of  the 


Table  1 

Names  and  locations  of  high  altitude  passes  on  the  border  between  Switzerland  and 
Italy  which  were  glacierized  in  1973. 


Number 

Name 

Latitude  (N) 

Longitude  (E) 

Altitude  (m) 

1 

Petit  Col  Ferret 

45°  53'  58" 

7°  4'  9" 

2490 

2 

Col  d'Amiante 

45°  55'  6" 

7°  18'  9" 

3319 

3 

Col  de  la  Balme 

45°  54'  7" 

7°  22'  27" 

3321 

4 

Col  du  Petit  Mont  Collon 

45°  57’  42" 

7°  29'  11" 

3292 

5 

Col  Collon 

45°  57'  41" 

7°  30'  51" 

3087 

6 

Col  des  Bouquetins 

45°  59’  13" 

7°  33'  36" 

3357 

7 

Col  de  la  Tete  Blanche 

45°  59'  33" 

7°  34'  53" 

3579 

8 

Tiefmattenjoch 

45°  58’  22" 

7°  35'  13" 

3543 

9 

Breuiljoch 

45°  58'  18" 

7°  40'  18" 

3313 

10 

Theodulpass 

45°  56'  38" 

7°  42’  35" 

3301 

11 

Passo  di  Ventina  Nord 

45°  56’  4" 

7°  42'  39" 

3450 

12 

Breithornpass  (south) 

46°  14'  36" 

8°  5’  12" 

3368 

13 

Zwillingsjoch 

45°  55’  36" 

7°  47'  24" 

3845 

14 

Felikjoch 

45°  55'  5" 

7°  48'  11" 

4066 

15 

Lisjoch 

45°  55'  18" 

7°  51'  12" 

4169 

16 

Neues  Weisstor 

45°  59'  20" 

7°  54’  4" 

3509 

17 

Seewjinenliicke 

45°  59’  53" 

7°  57'  22" 

3095 

18 

Tossenjoch 

46°  7'  31" 

8°  3’  22" 

2923 

19 

Breithornpass  (east) 

45°  55’  56" 

7°  44'  34" 

3845 
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Fig.  2.  Visualization  of  the  locational  analysis  processing  steps:  (a)  calculation  of  least  cost  paths  with  pass  number  corresponding  to  Table  1,  (b)  calculation  of  buffers  around  each 
path,  (c)  weighted  slope  values  derived  from  DEM,  (d)  weighted  glacier  thickness  layer,  (e)  selected  close-up  of  slope  multiplied  by  thickness,  and  (f)  the  final  weighted  locational 
analysis  layer  with  pass  numbers  (see  Table  1 ). 


terrain,  and  where  archaeological  remains  might  be  located  based 
on  glacial  characteristics  and  ice  thickness.  The  input  layers  are 
described  in  the  following  paragraphs. 

Buffers  were  constructed  around  each  path  in  100,  250,  500, 
750,  and  1000  m  intervals  on  each  side  of  the  path  to  represent  the 
notion  that  archaeological  potential  decreases  with  distance  from 
the  paths.  The  zones  within  the  100  m  buffers  were  assumed  to 
have  the  highest  archaeological  potential  and  weighted  with  the 
value  of  5,  while  the  zones  located  within  the  1  km  buffer  were 
assumed  to  have  to  lowest  potential  and  weighted  with  the  value  of 
1  (Table  2).  The  values  in  between  are  listed  in  Table  2. 

The  1000  m  buffers  calculated  above  were  used  to  define  the 
study  area  surrounding  the  paths  from  which  to  calculate  slope 
values.  The  extract  by  mask  tool  was  used  to  isolate  the  study  area 


for  the  DEM  and  the  slope  tool  was  used  to  calculate  the  steepness 
of  the  terrain  (Fig.  2c).  The  slope  values  were  calculated  and 
weighted  based  on  the  potential  of  finding  archaeological  remains. 
Slopes  greater  than  40°  were  given  a  value  of  1,  thus  very  low 


Table  2 

Weight  values  for  the  layers  used  in  the  locational  analysis. 


Weights 

Distance  from  LCPs  (m) 

Slope  (°) 

Ice  thickness  (m) 

5 

0-100 

0-10 

0-25 

4 

100-250 

10-20 

25-50 

3 

250-500 

20-30 

50-75 

2 

500-750 

30-40 

75-100 

1 

750-1000 

>40 

>100 
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archaeological  potential,  as  they  are  difficult  to  climb  and  would 
probably  have  been  avoided.  Furthermore,  at  steep  slopes  archae¬ 
ological  remains  are  more  likely  to  be  washed  out  by  erosion. 
Slopes  between  0  and  10°  were  determined  to  be  the  easiest  to 
walk  across  and  thus  were  given  the  highest  potential  value  of  5. 
Slopes  between  10  and  40°  were  assigned  with  a  respective  linearly 
increasing  potential  value  for  10°  classes  (Table  2). 

Next,  the  glacier  ice  thickness  layer  was  weighted  from  low  to 
high  archaeological  potential.  The  ice  thickness  range  between 
0  and  25  m  was  assigned  a  weight  of  5  as  those  are  the  areas  with 
the  highest  potential  now  and  in  future  years  because  they  will 
likely  be  the  first  to  become  ice-free  (Table  2).  Even  more  impor¬ 
tantly,  archeological  remains  are  only  likely  to  be  preserved  below 
relatively  thin  ice.  Thicker  ice  tends  to  flow  faster  and  is  generally 
more  destructive  in  terms  of  bedrock  erosion.  The  classes  with 
25—50  m,  50—75  m,  and  75—100  m  were  weighted  with  4,  3,  and  2, 
respectively  (Fig.  2d). 

The  weighted  layers  for  distance  from  LCPs,  slope,  and  thickness 
were  multiplied  together  to  obtain  one  layer  containing  all  possible 
value  combinations.  The  results  ranged  from  2  to  125  in  29  different 
classes  (Fig.  2e).  As  a  final  step  in  the  locational  analysis,  these  classes 
were  combined  into  five  final  potential  categories  using  the  Natural 
Breaks  classification  scheme  which  arranges  classes  into  “natural” 
objectively  selected  categories  (Fig.  2f)  (Jenks  and  Caspall,  1971). 

4.3.  Glaciological  modelling 

Future  changes  in  glacier  coverage  over  the  entire  Pennine  Alps 
were  assessed  by  a  combination  of  different  glaciological  models  at 
high  spatial  resolution.  First,  the  current  glacier  ice  thickness  dis¬ 
tribution  was  derived  using  the  glacier  outlines  from  Fischer  et  al. 
(2014)  for  Switzerland  and  from  Paul  et  al.  (2011)  for  Italy  based 
on  the  approach  by  Huss  and  Farinotti  (2012).  Next,  surface  mass 
balance  and  3D  glacier  geometry  change  were  modelled  transiently 
for  50  Swiss  glaciers  from  2010  to  2100  based  on  a  detailed  glacier 
model  (Huss  et  al.,  2010a).  The  model  runs  at  daily  resolution  on  a 
25  m  grid  and  takes  into  account  snow  accumulation  distribution, 
the  influence  of  radiation  on  ice  melting,  and  calculated  glacier 
retreat  based  on  a  mass-conservation  approach.  Model  calibration 
and  validation  for  the  50  investigated  glaciers  was  achieved  with  a 
variety  of  field  data  covering  the  entire  20th  century  (Huss  et  al., 
2010b).  For  calculating  future  glacier  change,  we  chose  to  use  one 
single  regional  climate  scenario  for  simplicity  although  the  pro¬ 
jected  evolution  of  meteorological  variables  is  subject  considerable 
uncertainties.  Seasonal  changes  in  air  temperature  and  precipita¬ 
tion  as  projected  by  the  Eidgenossische  Technische  Hochschule 
Zurich  (ETHZ,  Swiss  Federal  Institute  of  Technology)  RCM  were 
used  as  inputs  into  the  model.  This  climate  model  was  driven  by  the 
Al  B  C02-emission  scenario  (Nakicenovic,  2000).  Until  2100,  a  mean 
annual  air  temperature  rise  of  +4.7  °C  relative  to  1980—2009  is 
expected  for  the  study  region  and  precipitation  is  found  to  increase 
in  winter  but  to  decrease  in  summer  (CH2014-Impacts,  2014). 
Finally,  we  extrapolated  annual  mass  balance  from  the  50  glaciers 
to  every  glacier  in  the  Pennine  Alps  (Huss,  2012).  Thus,  for  each 
glacier,  a  glacier-specific  transient  annual  series  of  the  glacier  mass 
budget  was  obtained  which  was  used  to  drive  the  glacier  retreat 
model  (Huss  et  al.,  2010a).  From  the  transient  model  runs,  we 
extracted  glacier  ice  coverage  for  10-year  time  steps  between  2020 
and  2100.  These  glacier  masks  were  overlaid  onto  the  results  of  the 
locational  analysis. 

4.4.  GlaciArch 

In  this  step,  past  (1850  and  1973),  current  (2010),  and  selected 
future  (2030,  2060,  and  2090)  glacier  extents  were  overlaid  onto 


locational  analysis  results  to  create  the  GlaciArch  predictive  model. 
The  current  archaeological  potential  of  a  region  was  assessed 
differently  than  that  of  the  future  potential.  Current  archaeological 
potential  is  considered  to  exist  in  the  regions  that  have  been 
deglacierized  since  1973;  that  is,  between  the  1973  and  2010  ex¬ 
tents.  Those  are  the  general  areas  where  archaeological  remains 
could  be  currently  located  based  on  the  principles  of  glacier  dy¬ 
namics.  Future  archaeological  potential  is  ultimately  gauged  by 
comparing  the  modelled  glacier  extents  for  2030,  2060,  and  2090, 
to  the  results  obtained  by  the  locational  analysis. 

5.  Results  and  discussion 

5.3.  Locational  analysis 

The  locational  analysis  results  defined  regions  which  were  gla- 
cierized,  or  recently  deglacierized,  located  near  LCPs,  in  areas  with 
less  than  40°  slope,  and  where  ice  thickness  is  at  a  minimum.  The 
consideration  of  currently  glacierized  or  recently  deglacierized 
passes  is  important  when  dealing  with  glacial  archaeological  re¬ 
mains,  as  previously  suggested.  The  calculation  of  LCPs  enabled  a 
better  understanding  about  how  people  might  have  travelled  from 
one  location  to  another  based  on  the  principles  of  walking  across 
different  landscapes.  Although  prehistoric  and  historic  landscapes 
are  often  uncertain,  paleoecological  research  gives  a  good  general 
understanding  about  past  landscapes  (Berthel  et  al.,  2012;  Tinner 
and  Theurillat,  2003;  Wick  and  Tinner,  1997).  The  final  product  of 
the  locational  analysis  displays  the  overlapping  areas  of  the 
weighted  distance  from  LCPs,  slope,  and  ice  thickness  layers  in  the 
range  of  5  (high  potential)  to  1  (low  potential)  for  the  region 
(Fig.  3a).  The  results  of  locational  analysis  alone  reduced  a  region  of 
over  4500  km2  to  a  114  km2  area  of  interest,  8.16  km2  of  which  is 
considered  as  high  potential  (Table  3).  Similar  to  the  work  con¬ 
ducted  by  Dixon  et  al.  (2005)  and  Andrews  et  al.  (2012),  locational 
analysis  provided  a  means  to  define  small,  manageable  regions  for 
glacial  archaeological  prospection.  The  areas  are  more  finely 
delimited  in  the  final  GlaciArch  model  (Section  5.3). 

5.2.  Glaciological  modelling 

Future  glacier  extents  for  10-year  increments  between  2020  and 
2100  in  Switzerland  and  in  Italy  were  calculated  (shown  for  2030, 
2060,  and  2090  in  Fig.  3b).  This  regional  scale  modelling  method 
provided  a  high  resolution  projection  of  future  glacier  extents.  The 
total  area  of  glaciers  in  2010  in  the  Pennine  Alps  was  446  km2  and 
calculated  to  decrease  in  future  years.  For  example,  a  reduction  of 
37%— 280  km2  in  2030,  80%  to  91  km2  in  2060,  and  93%  to  30  km2  in 
2090  was  modelled  based  on  the  climate  scenario  used.  In  this 
study  we  did  not  assess  the  impact  of  different  assumptions  on 
future  climate  evolution  and  other  glaciological  model  un¬ 
certainties  on  the  results.  However,  Addor  et  al.  (in  press)  showed 
that  the  C02-emission  storyline  until  2100  only  had  a  small  effect 
on  calculated  total  glacier  area. 

5.3.  GlaciArch 

The  results  of  the  GlaciArch  model  show  current  and  future 
areas  of  archaeological  potential  spread  over  the  Pennine  Alps  re¬ 
gion  (see  supplementary  maps  covering  the  entire  region).  As 
mentioned  in  Section  5.2,  the  locational  analysis  results  provided  a 
broadly  defined  research  area.  The  addition  of  glaciological 
modelling  results  allows  the  delineations  to  be  further  defined.  For 
example,  between  2010  and  2030,  the  total  area  of  interest  is 
30.1  km2  and  the  high  potential  region  is  3.24  km2,  compared  to  the 
decreased  areas  between  2060  and  2090  with  13.7  km2  and 
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Fig.  3.  Results  from  (a)  locational  analysis  and  (b)  glaciological  modelling  from  the  whole  study  area.  The  Swiss  Glacier  Inventory  for  2010  is  shown  along  with  the  modelled  glacier 
extents  for  2010  on  the  Italian  side  of  the  border.  The  projected  glacier  extents  for  2030, 2060,  and  2090  are  shown  for  both  sides  of  the  border  of  the  Pennine  Alps.  The  boundaries 
for  Figs.  1  and  2  are  shown  in  red.  Boundaries  SI  -S5  refer  to  additional  maps  not  discussed  in  the  text  which  can  be  found  in  the  supplementary  data  section.  (For  interpretation  of 
the  references  to  colour  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


0.58  km2,  for  area  of  interest  and  high  potential,  respectively  (see 
Table  3  for  all  values). 

Taking  a  closer  look  at  the  results,  the  area  surrounding  the 
Mont  Collon  (Fig.  4)  is  an  area  of  interest  due  to  its  archaeological 
significance  in  the  surrounding  regions,  although  no  glacial 
archaeological  finds  have  been  retrieved  from  that  site  to  date 
(Fig.  4).  In  1948,  a  Neolithic  tool  made  of  flint  was  located  on  the 
Plans  de  Bertol,  which  is  located  on  the  way  to  the  Col  Collon,  which 
leads  archaeologists  to  believe  that  it  was  a  pass  which  was  used  by 
humans  for  thousands  of  years  (Bezinge  and  Curdy,  1994,  1995; 
Sauter,  1950)  (Fig.  4).  Currently,  the  model  indicates  that  there  is 


high  glacial  archaeological  potential  on  the  margins  of  the  Haut 
Glacier  d'Arolla  and  the  Glacier  du  Mont  Collon,  as  well  as  the  area 
between  the  Petit  Mont  Collon  and  the  northern  section  of  the 
Glacier  d'Otemma,  just  south  of  the  Col  de  Charmontagne  (Fig.  4). 
For  the  future,  the  results  from  GlaciArch  show  that  between  2010 
and  2030,  there  will  be  high  potential  areas  on  the  margins  of  the 
tongue  of  the  Glacier  du  Mont  Collon,  as  well  as  the  Haut  Glacier 
d'Arolla.  Between  2030  and  2060,  the  Col  Collon  (5)  is  expected  to 
become  an  area  of  high  archaeological  potential  as  well  as  some 
regions  surrounding  the  Petit  Mont  Collon.  Sections  in  the  middle 
of  the  Glacier  du  Mont  Collon  and  to  the  north  of  the  southern 
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Table  3 

Locational  analysis  areas  calculated  for  2010  and  the  time  periods  between 
2010—2030,  2030-2060,  and  2060—2090.  The  high  glacial  archaeological  potential 
(value  5)  areas  are  also  calculated  for  each  time  period. 


Year(s) 

Total  area  (km2) 

High  potential  (km2) 

2010 

114 

8.16 

2010-2030 

30.1 

3.24 

2030-2060 

40.5 

2.31 

2060-2090 

13.7 

0.58 

section  of  the  Glacier  d'Otemma  also  show  high  potential.  Ac¬ 
cording  to  the  results,  it  appears  that  by  2060,  the  Haut  Glacier 
d'Arolla  will  have  almost  completely  disappeared.  By  2090,  very 
little  ice  will  remain,  however,  the  section  to  the  north  of  the  Col  du 
Petit  Mont  Collon  (4)  could  be  of  interest  at  that  time. 

Already  briefly  discussed  in  Section  3,  the  area  surrounding  the 
Theodulpass  will  now  be  revisited.  The  high  altitude  passes  of  in¬ 
terest  near  the  Theodulpass  resulting  from  this  analysis  are  the 
Breuiljoch  (9),  Theodulpass  (10),  Passo  di  Ventina  Nord  ( 11 ),  and  the 
Breithornpass  (12)  (Fig.  5).  From  2010  to  2030  there  are  regions  of 
high  archaeological  potential  on  the  extents  of  the  Furgg,  Oberer 
Theodul,  and  Valtournanche  glaciers,  as  well  as  on  the  Theodulpass 
and  Passo  de  Ventina  Nord.  Between  2030  and  2060,  the  east  side  of 
the  Oberer  Theodul  glacier  becomes  a  predominant  area  of  interest, 
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Fig.  4.  GlaciArch  results  around  the  Mont  Collon  which  includes  the  Col  de  Petit  Mont 
Collon  (4)  and  the  Col  Collon  (5).  The  Glacier  du  Mont  Collon,  Glacier  d'Otemma,  and 
Haut  glacier  d'Arolla  are  labelled  as  well  as  locations  mentioned  in  the  text:  Petit  Mont 
Collon,  Col  de  Charmontagne,  and  Plans  de  Bertol.  LCP  refers  to  least  cost  path. 


while  the  Furgg  glacier  has  decreased  archaeological  potential.  The 
area  south  of  the  Passo  de  Ventina  also  shows  high  potential  during 
that  time  period.  From  2060  to  2090  The  Furgg  and  Oberer  Theodul 
glaciers  are  predicted  to  almost  completely  disappear  therefore 
their  glacial  archaeological  potential  decreases  significantly. 

The  performance  of  this  model  has  yet  to  be  tested  in  the  field, 
however  the  results  seem  to  correspond  well  to  glaciological 
principles  and  the  few  glacial  archaeological  finds  already  located 
in  the  region,  for  example  at  the  Oberer  Theodul  site  (Figs.  4  and  5). 
In  theory,  high  archaeological  potential  is  expected  near  the  glacier 
margins  as  those  are  the  areas  with  the  thinnest  ice,  while  areas  of 
low  potential  should  occur  mid-glacier  where  the  ice  is  thickest.  An 
example  of  this  can  be  seen  in  Fig.  4  at  the  Haut  Glacier  d'Arolla; 
high  potential  exists  on  the  extents  of  the  glacier  tongue  and  po¬ 
tential  decreases  as  inward  movement  onto  the  glacier  continues. 
One  problem  with  this  model  is  that  it  is  difficult  to  convey 
inherently  dynamic  movement  on  a  static  map.  In  fact,  potential 
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Fig.  5.  GlaciArch  results  around  the  Theodulhorn.  Passes  listed  are  the  Breuiljoch  (9), 
Theodulpass  (10),  and  the  Passo  di  Ventina  Nord  (11).  The  Furgg  glacier,  and  Oberer 
and  Unterer  Theodul  glaciers  are  labelled.  LCP  refers  to  least  cost  path. 
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maps  should  be  calculated  for  each  year  to  obtain  a  greater  un¬ 
derstanding  about  the  true  and  temporally  varying  archaeological 
potential  of  the  region,  however  2D  mapping  contraints  do  not 
permit  dynamic  visualisation  of  this  type  of  results.  An  interactive 
interface  would  be  the  best  way  to  visualize  the  end  results. 

6.  Conclusion 

In  this  paper,  the  new  integrated  model  GlaciArch  was  used  to 
identify  areas  of  current  and  future  archaeological  interest  and 
potential  in  the  Pennine  Alps.  The  model  highlights  areas  which 
correspond  to  the  retrieval  of  archaeological  remains  based  on 
glaciological  principles  and  the  topographic  properties  of  the 
terrain.  Thus,  the  GlaciArch  results  can  be  used  as  a  decision  sup¬ 
port  tool  for  the  selection  of  glacial  archaeological  prospection 
sites.  The  definition  of  small  regions  of  high  glacial  archaeological 
potential  means  less  time,  effort,  and  money  spent  in  the  field  or  on 
flight  reconnaissance  missions.  By  combining  archaeological  and 
glaciological  methods  for  the  first  time,  a  new  perspective  has  been 
given  to  the  field  of  glacial  archaeology.  The  integration  of  loca¬ 
tional  analysis  and  regional  scale  glacier  modelling  proved  to  be 
beneficial  for  narrowing  down  a  large,  often  inaccessible,  and 
remote  study  region  to  identify  zones  of  archaeological  interest 
based  on  glaciological  characteristics  and  human  accessibility.  The 
glacier  modelling  results  forecast  a  93%  loss  in  area  by  2090.  With 
these  alarming  melting  rates,  immediate  focus  should  be  given  to 
high  archaeological  potential  areas  in  hopes  to  locate  and  recover 
possibly  irreplaceable,  culturally  significant  items.  In  order  to  pro¬ 
tect  and  conserve  these  exceptional  and  rare  relics,  further  multi¬ 
disciplinary  predictive  methods  should  be  developed  and 
employed  in  sensitive  areas  such  as  the  Pennine  Alps. 
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